Metropolis algorithm

ggplot2 is used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(tidyr)
library(gganimate)
library(ggforce)
library(MASS)
library(rprojroot)
library(rstan)
root<-has_file(".BDA_R_demos_root")$make_fix_file()

Parameters of a normal distribution used as a toy target distribution

y1 <- 0
y2 <- 0
r <- 0.8
S <- diag(2)
S[1, 2] <- r
S[2, 1] <- r

Metropolis proposal distribution scale

sp <- 0.3

Sample from the toy distribution to visualize 90% HPD interval with ggplot's stat_ellipse()

dft <- data.frame(mvrnorm(100000, c(0, 0), S))

see BDA3 p. 85 for how to compute HPD for multivariate normal in 2d-case contour for 90% HPD is an ellipse, whose semimajor axes can be computed from the eigenvalues of the covariance matrix scaled by a value selected to get ellipse match the density at the edge of 90% HPD. Angle of the ellipse could be computed from the eigenvectors, but since the marginals are same we know that angle is pi/4 Starting value of the chain

t1 <- -2.5
t2 <- 2.5

Number of iterations.

M <- 5000

Insert your own Metropolis sampling here

# Allocate memory for the sample
tt <- matrix(rep(0, 2*M), ncol = 2)
tt[1,] <- c(t1, t2)    # Save starting point
# For demonstration load pre-computed values
# Replace this with your algorithm!
# tt is a M x 2 array, with M draws of both theta_1 and theta_2
load(root("demos_ch11","demo11_2a.RData"))

The rest is for illustration Take the first 200 draws to illustrate how the sampler works

df100 <- data.frame(id=rep(1,100),
                    iter=1:100, 
                    th1 = tt[1:100, 1],
                    th2 = tt[1:100, 2],
                    th1l = c(tt[1, 1], tt[1:(100-1), 1]),
                    th2l = c(tt[1, 2], tt[1:(100-1), 2]))

Take the first 5000 observations after warmup of 50

s <- 5000
warm <- 500
dfs <- data.frame(th1 = tt[(warm+1):s, 1], th2 = tt[(warm+1):s, 2])

Remove warm-up period of 50 first draws later

# labels and frame indices for the plot
labs1 <- c('Draws', 'Steps of the sampler', '90% HPD')
p1 <- ggplot() +
  geom_jitter(data = df100, width=0.05, height=0.05,
              aes(th1, th2, group=id, color ='1'), alpha=0.3) +
  geom_segment(data = df100, aes(x = th1, xend = th1l, color = '2',
                                 y = th2, yend = th2l)) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '3'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('red', 'forestgreen','blue'), labels = labs1) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA, NA), linetype = c(0, 1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

The following generates a gif animation of the steps of the sampler (might take 10 seconds).

animate(p1 +   
          transition_reveal(along=iter) + 
          shadow_trail(0.01))
## Rendering [--------------------------------------------] at 4.5 fps ~ eta: 22s
## Rendering [>---------------------------------------------] at 5 fps ~ eta: 20s
## Rendering [>---------------------------------------------] at 5 fps ~ eta: 19s
## Rendering [=>------------------------------------------] at 5.1 fps ~ eta: 19s
## Rendering [=>------------------------------------------] at 5.2 fps ~ eta: 18s
## Rendering [==>-----------------------------------------] at 5.2 fps ~ eta: 18s
## Rendering [===>----------------------------------------] at 5.2 fps ~ eta: 18s
## Rendering [===>----------------------------------------] at 5.2 fps ~ eta: 17s
## Rendering [====>---------------------------------------] at 5.2 fps ~ eta: 17s
## Rendering [=====>--------------------------------------] at 5.3 fps ~ eta: 17s
## Rendering [=====>--------------------------------------] at 5.3 fps ~ eta: 16s
## Rendering [======>-------------------------------------] at 5.2 fps ~ eta: 16s
## Rendering [======>-------------------------------------] at 5.3 fps ~ eta: 16s
## Rendering [=======>------------------------------------] at 5.3 fps ~ eta: 16s
## Rendering [=======>------------------------------------] at 5.3 fps ~ eta: 15s
## Rendering [========>-----------------------------------] at 5.1 fps ~ eta: 16s
## Rendering [========>-----------------------------------] at 5.1 fps ~ eta: 15s
## Rendering [=========>----------------------------------] at 5.1 fps ~ eta: 15s
## Rendering [==========>---------------------------------] at 5.1 fps ~ eta: 15s
## Rendering [==========>---------------------------------] at 5.1 fps ~ eta: 14s
## Rendering [===========>--------------------------------] at 5.1 fps ~ eta: 14s
## Rendering [===========>--------------------------------] at 5.2 fps ~ eta: 14s
## Rendering [============>-------------------------------] at 5.2 fps ~ eta: 14s
## Rendering [=============>------------------------------] at 5.2 fps ~ eta: 13s
## Rendering [=============>------------------------------] at 5.1 fps ~ eta: 13s
## Rendering [==============>-----------------------------] at 5.1 fps ~ eta: 13s
## Rendering [===============>------------------------------] at 5 fps ~ eta: 13s
## Rendering [===============>----------------------------] at 5.1 fps ~ eta: 13s
## Rendering [===============>----------------------------] at 5.1 fps ~ eta: 12s
## Rendering [================>---------------------------] at 5.1 fps ~ eta: 12s
## Rendering [=================>--------------------------] at 5.1 fps ~ eta: 12s
## Rendering [=================>--------------------------] at 5.1 fps ~ eta: 11s
## Rendering [==================>-------------------------] at 5.1 fps ~ eta: 11s
## Rendering [===================>------------------------] at 5.1 fps ~ eta: 11s
## Rendering [=====================>------------------------] at 5 fps ~ eta: 11s
## Rendering [=====================>------------------------] at 5 fps ~ eta: 10s
## Rendering [=====================>----------------------] at 4.9 fps ~ eta: 10s
## Rendering [=====================>----------------------] at 4.8 fps ~ eta: 10s
## Rendering [======================>---------------------] at 4.8 fps ~ eta: 10s
## Rendering [======================>---------------------] at 4.7 fps ~ eta: 10s
## Rendering [=======================>--------------------] at 4.7 fps ~ eta: 10s
## Rendering [========================>-------------------] at 4.7 fps ~ eta: 9s
## Rendering [=========================>------------------] at 4.7 fps ~ eta: 9s
## Rendering [=========================>------------------] at 4.6 fps ~ eta: 9s
## Rendering [==========================>-----------------] at 4.6 fps ~ eta: 8s
## Rendering [==========================>-----------------] at 4.7 fps ~ eta: 8s
## Rendering [===========================>----------------] at 4.7 fps ~ eta: 8s
## Rendering [============================>---------------] at 4.7 fps ~ eta: 8s
## Rendering [============================>---------------] at 4.6 fps ~ eta: 7s
## Rendering [=============================>--------------] at 4.6 fps ~ eta: 7s
## Rendering [==============================>-------------] at 4.6 fps ~ eta: 6s
## Rendering [===============================>------------] at 4.6 fps ~ eta: 6s
## Rendering [================================>-----------] at 4.6 fps ~ eta: 6s
## Rendering [================================>-----------] at 4.5 fps ~ eta: 6s
## Rendering [================================>-----------] at 4.5 fps ~ eta: 5s
## Rendering [=================================>----------] at 4.4 fps ~ eta: 5s
## Rendering [==================================>---------] at 4.4 fps ~ eta: 5s
## Rendering [===================================>--------] at 4.3 fps ~ eta: 4s
## Rendering [====================================>-------] at 4.3 fps ~ eta: 4s
## Rendering [====================================>-------] at 4.3 fps ~ eta: 3s
## Rendering [=====================================>------] at 4.3 fps ~ eta: 3s
## Rendering [======================================>-----] at 4.4 fps ~ eta: 3s
## Rendering [=======================================>----] at 4.4 fps ~ eta: 2s
## Rendering [========================================>---] at 4.4 fps ~ eta: 2s
## Rendering [========================================>---] at 4.4 fps ~ eta: 1s
## Rendering [=========================================>--] at 4.4 fps ~ eta: 1s
## Rendering [==========================================>-] at 4.4 fps ~ eta: 1s
## Rendering [==========================================>-] at 4.4 fps ~ eta: 0s
## Rendering [===========================================>] at 4.4 fps ~ eta: 0s
## Rendering [============================================] at 4.4 fps ~ eta: 0s

Plot the final frame

p1

show 1000 draws after the warm-up

labs2 <- c('Draws', '90% HPD')
ggplot() +
  geom_point(data = dfs[1:1000,],
             aes(th1, th2, color = '1'), alpha = 0.3) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '2'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('steelblue', 'blue'), labels = labs2) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA), linetype = c(0, 1), alpha = c(1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

show 4500 draws after the warm-up

labs2 <- c('Draws', '90% HPD')
ggplot() +
  geom_point(data = dfs,
             aes(th1, th2, color = '1'), alpha = 0.3) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '2'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('steelblue', 'blue'), labels = labs2) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA), linetype = c(0, 1), alpha = c(1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Convergence diagnostics

samp <- tt
dim(samp) <- c(dim(tt),1)
samp <- aperm(samp, c(1, 3, 2))
res<-monitor(samp, probs = c(0.25, 0.5, 0.75), digits_summary = 2)
## Inference for the input samples (1 chains: each with iter = 5000; warmup = 2500):
## 
##      Q5 Q50 Q95 Mean SD  Rhat Bulk_ESS Tail_ESS
## V1 -1.8   0 1.6 -0.1  1     1      101      152
## V2 -1.7   0 1.7  0.0  1     1       92      105
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).
neff <- res[,'n_eff']
# both theta have owen neff, but for plotting these are so close to each
# other, so that single relative efficiency value is used
reff <- mean(neff/(s/2))

Visual convergence diagnostics

Collapse the data frame with row numbers augmented into key-value pairs for visualizing the chains

dfb <- dfs
sb <- s-warm
dfch <- within(dfb, iter <- 1:sb) %>% gather(grp, value, -iter)

Another data frame for visualizing the estimate of the autocorrelation function

nlags <- 50
dfa <- sapply(dfb, function(x) acf(x, lag.max = nlags, plot = F)$acf) %>%
  data.frame(iter = 0:(nlags)) %>% gather(grp, value, -iter)

A third data frame to visualize the cumulative averages and the 95% intervals

dfca <- (cumsum(dfb) / (1:sb)) %>%
  within({iter <- 1:sb
  uppi <-  1.96/sqrt(1:sb)
  upp <- 1.96/(sqrt(1:sb*reff))}) %>%
  gather(grp, value, -iter)

Visualize the chains

ggplot(data = dfch) +
  geom_line(aes(iter, value, color = grp)) +
  labs(title = 'Trends') +
  scale_color_discrete(labels = c('theta1','theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the autocorrelation function

ggplot(data = dfa) +
  geom_line(aes(iter, value, color = grp)) +
  geom_hline(aes(yintercept = 0)) +
  labs(title = 'Autocorrelation function') +
  scale_color_discrete(labels = c('theta1', 'theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the Monte Carlo error estimates

# labels
labs3 <- c('theta1', 'theta2',
           '95% interval for MCMC error',
           '95% interval for independent MC')
ggplot() +
  geom_line(data = dfca, aes(iter, value, color = grp, linetype = grp)) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb*reff)), linetype = 2) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb)), linetype = 3) +
  geom_hline(aes(yintercept = 0)) +
  coord_cartesian(ylim = c(-1.5, 1.5), xlim = c(0,4000)) +
  labs(title = 'Cumulative averages') +
  scale_color_manual(values = c('red','blue',rep('black', 2)), labels = labs3) +
  scale_linetype_manual(values = c(1, 1, 2, 3), labels = labs3) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Same again with r=0.99 Parameters of a normal distribution used as a toy target distribution

y1 <- 0
y2 <- 0
r <- 0.99
S <- diag(2)
S[1, 2] <- r
S[2, 1] <- r

Metropolis proposal distribution scale

sp <- 0.3

Sample from the toy distribution to visualize 90% HPD interval with ggplot's stat_ellipse()

dft <- data.frame(mvrnorm(100000, c(0, 0), S))

see BDA3 p. 85 for how to compute HPD for multivariate normal in 2d-case contour for 90% HPD is an ellipse, whose semimajor axes can be computed from the eigenvalues of the covariance matrix scaled by a value selected to get ellipse match the density at the edge of 90% HPD. Angle of the ellipse could be computed from the eigenvectors, but since the marginals are same we know that angle is pi/4 Starting value of the chain

t1 <- -2.5
t2 <- 2.5

Number of iterations.

M <- 5000

Insert your own Metropolis sampling here

# Allocate memory for the sample
tt <- matrix(rep(0, 2*M), ncol = 2)
tt[1,] <- c(t1, t2)    # Save starting point
# For demonstration load pre-computed values
# Replace this with your algorithm!
# tt is a M x 2 array, with M draws of both theta_1 and theta_2
load(root("demos_ch11","demo11_2b.RData"))

The rest is for illustration Take the first 200 draws to illustrate how the sampler works

df100 <- data.frame(id=rep(1,100),
                    iter=1:100, 
                    th1 = tt[1:100, 1],
                    th2 = tt[1:100, 2],
                    th1l = c(tt[1, 1], tt[1:(100-1), 1]),
                    th2l = c(tt[1, 2], tt[1:(100-1), 2]))

Take the first 5000 observations after warmup of 50

s <- 5000
warm <- 500
dfs <- data.frame(th1 = tt[(warm+1):s, 1], th2 = tt[(warm+1):s, 2])

Remove warm-up period of 50 first draws later

# labels and frame indices for the plot
labs1 <- c('Draws', 'Steps of the sampler', '90% HPD')
p1 <- ggplot() +
  geom_jitter(data = df100, width=0.05, height=0.05,
             aes(th1, th2, group=id, color ='1'), alpha=0.3) +
  geom_segment(data = df100, aes(x = th1, xend = th1l, color = '2',
                                 y = th2, yend = th2l)) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '3'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('red', 'forestgreen','blue'), labels = labs1) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA, NA), linetype = c(0, 1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

The following generates a gif animation of the steps of the sampler (might take 10 seconds).

animate(p1 +   
          transition_reveal(along=iter) + 
        shadow_trail(0.01))
## Rendering [>-------------------------------------------] at 5.6 fps ~ eta: 18s
## Rendering [>-------------------------------------------] at 4.5 fps ~ eta: 21s
## Rendering [=>--------------------------------------------] at 4 fps ~ eta: 24s
## Rendering [=>------------------------------------------] at 3.6 fps ~ eta: 26s
## Rendering [==>-----------------------------------------] at 3.2 fps ~ eta: 29s
## Rendering [==>-----------------------------------------] at 3.3 fps ~ eta: 28s
## Rendering [===>----------------------------------------] at 3.5 fps ~ eta: 27s
## Rendering [===>----------------------------------------] at 3.4 fps ~ eta: 27s
## Rendering [====>---------------------------------------] at 3.3 fps ~ eta: 27s
## Rendering [====>---------------------------------------] at 3.3 fps ~ eta: 26s
## Rendering [=====>--------------------------------------] at 3.4 fps ~ eta: 26s
## Rendering [=====>--------------------------------------] at 3.5 fps ~ eta: 25s
## Rendering [======>-------------------------------------] at 3.5 fps ~ eta: 24s
## Rendering [======>-------------------------------------] at 3.6 fps ~ eta: 23s
## Rendering [======>-------------------------------------] at 3.7 fps ~ eta: 22s
## Rendering [=======>------------------------------------] at 3.8 fps ~ eta: 22s
## Rendering [=======>------------------------------------] at 3.8 fps ~ eta: 21s
## Rendering [========>-----------------------------------] at 3.8 fps ~ eta: 21s
## Rendering [========>-----------------------------------] at 3.9 fps ~ eta: 20s
## Rendering [=========>----------------------------------] at 3.9 fps ~ eta: 20s
## Rendering [==========>-----------------------------------] at 4 fps ~ eta: 19s
## Rendering [==========>---------------------------------] at 4.1 fps ~ eta: 19s
## Rendering [==========>---------------------------------] at 4.1 fps ~ eta: 18s
## Rendering [===========>--------------------------------] at 4.1 fps ~ eta: 18s
## Rendering [===========>--------------------------------] at 4.2 fps ~ eta: 17s
## Rendering [============>-------------------------------] at 4.2 fps ~ eta: 17s
## Rendering [=============>------------------------------] at 4.2 fps ~ eta: 16s
## Rendering [=============>------------------------------] at 4.3 fps ~ eta: 16s
## Rendering [==============>-----------------------------] at 4.3 fps ~ eta: 16s
## Rendering [==============>-----------------------------] at 4.3 fps ~ eta: 15s
## Rendering [===============>----------------------------] at 4.3 fps ~ eta: 15s
## Rendering [===============>----------------------------] at 4.3 fps ~ eta: 14s
## Rendering [================>---------------------------] at 4.4 fps ~ eta: 14s
## Rendering [=================>--------------------------] at 4.4 fps ~ eta: 14s
## Rendering [=================>--------------------------] at 4.4 fps ~ eta: 13s
## Rendering [=================>--------------------------] at 4.5 fps ~ eta: 13s
## Rendering [==================>-------------------------] at 4.5 fps ~ eta: 13s
## Rendering [==================>-------------------------] at 4.5 fps ~ eta: 12s
## Rendering [===================>------------------------] at 4.5 fps ~ eta: 12s
## Rendering [====================>-----------------------] at 4.5 fps ~ eta: 12s
## Rendering [====================>-----------------------] at 4.5 fps ~ eta: 11s
## Rendering [=====================>----------------------] at 4.5 fps ~ eta: 11s
## Rendering [=====================>----------------------] at 4.6 fps ~ eta: 11s
## Rendering [======================>---------------------] at 4.6 fps ~ eta: 11s
## Rendering [======================>---------------------] at 4.5 fps ~ eta: 10s
## Rendering [=======================>--------------------] at 4.5 fps ~ eta: 10s
## Rendering [========================>-------------------] at 4.5 fps ~ eta: 10s
## Rendering [=========================>------------------] at 4.5 fps ~ eta: 9s
## Rendering [==========================>-----------------] at 4.6 fps ~ eta: 9s
## Rendering [==========================>-----------------] at 4.6 fps ~ eta: 8s
## Rendering [===========================>----------------] at 4.6 fps ~ eta: 8s
## Rendering [===========================>----------------] at 4.5 fps ~ eta: 8s
## Rendering [============================>---------------] at 4.5 fps ~ eta: 8s
## Rendering [============================>---------------] at 4.5 fps ~ eta: 7s
## Rendering [=============================>--------------] at 4.6 fps ~ eta: 7s
## Rendering [==============================>-------------] at 4.6 fps ~ eta: 7s
## Rendering [==============================>-------------] at 4.6 fps ~ eta: 6s
## Rendering [===============================>------------] at 4.6 fps ~ eta: 6s
## Rendering [================================>-----------] at 4.6 fps ~ eta: 6s
## Rendering [================================>-----------] at 4.6 fps ~ eta: 5s
## Rendering [=================================>----------] at 4.6 fps ~ eta: 5s
## Rendering [=================================>----------] at 4.5 fps ~ eta: 5s
## Rendering [==================================>---------] at 4.5 fps ~ eta: 5s
## Rendering [==================================>---------] at 4.5 fps ~ eta: 4s
## Rendering [===================================>--------] at 4.5 fps ~ eta: 4s
## Rendering [====================================>-------] at 4.5 fps ~ eta: 4s
## Rendering [====================================>-------] at 4.5 fps ~ eta: 3s
## Rendering [=====================================>------] at 4.5 fps ~ eta: 3s
## Rendering [======================================>-----] at 4.5 fps ~ eta: 3s
## Rendering [======================================>-----] at 4.5 fps ~ eta: 2s
## Rendering [=======================================>----] at 4.5 fps ~ eta: 2s
## Rendering [========================================>---] at 4.5 fps ~ eta: 2s
## Rendering [========================================>---] at 4.6 fps ~ eta: 1s
## Rendering [=========================================>--] at 4.6 fps ~ eta: 1s
## Rendering [==========================================>-] at 4.6 fps ~ eta: 1s
## Rendering [==========================================>-] at 4.5 fps ~ eta: 0s
## Rendering [===========================================>] at 4.5 fps ~ eta: 0s
## Rendering [============================================] at 4.5 fps ~ eta: 0s

Plot the final frame

p1

show 1000 draws after the warm-up

labs2 <- c('Draws', '90% HPD')
ggplot() +
  geom_point(data = dfs[1:1000,],
             aes(th1, th2, color = '1'), alpha = 0.3) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '2'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('steelblue', 'blue'), labels = labs2) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA), linetype = c(0, 1), alpha = c(1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

show 4500 draws after the warm-up

labs2 <- c('Draws', '90% HPD')
ggplot() +
  geom_point(data = dfs,
             aes(th1, th2, color = '1'), alpha = 0.3) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '2'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('steelblue', 'blue'), labels = labs2) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA), linetype = c(0, 1), alpha = c(1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Convergence diagnostics

samp <- tt
dim(samp) <- c(dim(tt),1)
samp <- aperm(samp, c(1, 3, 2))
res<-monitor(samp, probs = c(0.25, 0.5, 0.75), digits_summary = 2)
## Inference for the input samples (1 chains: each with iter = 5000; warmup = 2500):
## 
##      Q5  Q50 Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## V1 -1.7 -0.4 1.5 -0.2 1.1  1.02       24       50
## V2 -1.6 -0.4 1.5 -0.2 1.1  1.03       24       51
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).
neff <- res[,'n_eff']
# both theta have owen neff, but for plotting these are so close to each
# other, so that single relative efficiency value is used
reff <- mean(neff/(s/2))

Visual convergence diagnostics

Collapse the data frame with row numbers augmented into key-value pairs for visualizing the chains

dfb <- dfs
sb <- s-warm
dfch <- within(dfb, iter <- 1:sb) %>% gather(grp, value, -iter)

Another data frame for visualizing the estimate of the autocorrelation function

nlags <- 100
dfa <- sapply(dfb, function(x) acf(x, lag.max = nlags, plot = F)$acf) %>%
  data.frame(iter = 0:(nlags)) %>% gather(grp, value, -iter)

A third data frame to visualize the cumulative averages and the 95% intervals

dfca <- (cumsum(dfb) / (1:sb)) %>%
  within({iter <- 1:sb
          uppi <-  1.96/sqrt(1:sb)
          upp <- 1.96/(sqrt(1:sb*reff))}) %>%
  gather(grp, value, -iter)

Visualize the chains

ggplot(data = dfch) +
  geom_line(aes(iter, value, color = grp)) +
  labs(title = 'Trends') +
  scale_color_discrete(labels = c('theta1','theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the autocorrelation function

ggplot(data = dfa) +
  geom_line(aes(iter, value, color = grp)) +
  geom_hline(aes(yintercept = 0)) +
  labs(title = 'Autocorrelation function') +
  scale_color_discrete(labels = c('theta1', 'theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the Monte Carlo error estimates

# labels
labs3 <- c('theta1', 'theta2',
           '95% interval for MCMC error',
           '95% interval for independent MC')
ggplot() +
  geom_line(data = dfca, aes(iter, value, color = grp, linetype = grp)) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb*reff)), linetype = 2) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb)), linetype = 3) +
  geom_hline(aes(yintercept = 0)) +
  coord_cartesian(ylim = c(-1.5, 1.5), xlim = c(0,4000)) +
  labs(title = 'Cumulative averages') +
  scale_color_manual(values = c('red','blue',rep('black', 2)), labels = labs3) +
  scale_linetype_manual(values = c(1, 1, 2, 3), labels = labs3) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Same again with sp = 1.5

sp = 1.5

Insert your own Metropolis sampling here

# Allocate memory for the sample
tt <- matrix(rep(0, 2*M), ncol = 2)
tt[1,] <- c(t1, t2)    # Save starting point
# For demonstration load pre-computed values
# Replace this with your algorithm!
# tt is a M x 2 array, with M draws of both theta_1 and theta_2
load(root("demos_ch11","demo11_2c.RData"))

The rest is for illustration Take the first 200 draws to illustrate how the sampler works

df100 <- data.frame(id=rep(1,100),
                    iter=1:100, 
                    th1 = tt[1:100, 1],
                    th2 = tt[1:100, 2],
                    th1l = c(tt[1, 1], tt[1:(100-1), 1]),
                    th2l = c(tt[1, 2], tt[1:(100-1), 2]))

Take the first 5000 observations after warmup of 50

s <- 5000
warm <- 500
dfs <- data.frame(th1 = tt[(warm+1):s, 1], th2 = tt[(warm+1):s, 2])

Remove warm-up period of 50 first draws later

# labels and frame indices for the plot
labs1 <- c('Draws', 'Steps of the sampler', '90% HPD')
p1 <- ggplot() +
  geom_jitter(data = df100, width=0.05, height=0.05,
             aes(th1, th2, group=id, color ='1'), alpha=0.3) +
  geom_segment(data = df100, aes(x = th1, xend = th1l, color = '2',
                                 y = th2, yend = th2l)) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '3'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('red', 'forestgreen','blue'), labels = labs1) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA, NA), linetype = c(0, 1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

The following generates a gif animation of the steps of the sampler (might take 10 seconds).

animate(p1 +   
          transition_reveal(along=iter) + 
        shadow_trail(0.01))
## Rendering [--------------------------------------------] at 3.3 fps ~ eta: 30s
## Rendering [>-------------------------------------------] at 2.8 fps ~ eta: 35s
## Rendering [>-------------------------------------------] at 2.9 fps ~ eta: 33s
## Rendering [=>------------------------------------------] at 2.8 fps ~ eta: 34s
## Rendering [==>-----------------------------------------] at 2.8 fps ~ eta: 33s
## Rendering [===>----------------------------------------] at 2.8 fps ~ eta: 33s
## Rendering [===>----------------------------------------] at 2.9 fps ~ eta: 32s
## Rendering [===>----------------------------------------] at 2.9 fps ~ eta: 31s
## Rendering [====>---------------------------------------] at 2.9 fps ~ eta: 31s
## Rendering [====>---------------------------------------] at 2.9 fps ~ eta: 30s
## Rendering [=====>--------------------------------------] at 2.9 fps ~ eta: 30s
## Rendering [======>-------------------------------------] at 2.9 fps ~ eta: 29s
## Rendering [=======>------------------------------------] at 2.9 fps ~ eta: 28s
## Rendering [========>-----------------------------------] at 2.9 fps ~ eta: 28s
## Rendering [========>-----------------------------------] at 2.9 fps ~ eta: 27s
## Rendering [=========>----------------------------------] at 2.9 fps ~ eta: 27s
## Rendering [=========>----------------------------------] at 2.8 fps ~ eta: 27s
## Rendering [==========>---------------------------------] at 2.8 fps ~ eta: 27s
## Rendering [===========>--------------------------------] at 2.8 fps ~ eta: 26s
## Rendering [============>-------------------------------] at 2.8 fps ~ eta: 26s
## Rendering [============>-------------------------------] at 2.8 fps ~ eta: 25s
## Rendering [=============>------------------------------] at 2.8 fps ~ eta: 25s
## Rendering [==============>-----------------------------] at 2.8 fps ~ eta: 24s
## Rendering [==============>-----------------------------] at 2.8 fps ~ eta: 23s
## Rendering [===============>----------------------------] at 2.9 fps ~ eta: 22s
## Rendering [================>---------------------------] at 2.9 fps ~ eta: 21s
## Rendering [=================>----------------------------] at 3 fps ~ eta: 20s
## Rendering [==================>---------------------------] at 3 fps ~ eta: 20s
## Rendering [==================>---------------------------] at 3 fps ~ eta: 19s
## Rendering [==================>-------------------------] at 3.1 fps ~ eta: 19s
## Rendering [==================>-------------------------] at 3.1 fps ~ eta: 18s
## Rendering [===================>------------------------] at 3.1 fps ~ eta: 18s
## Rendering [===================>------------------------] at 3.1 fps ~ eta: 17s
## Rendering [====================>-----------------------] at 3.1 fps ~ eta: 17s
## Rendering [====================>-----------------------] at 3.2 fps ~ eta: 16s
## Rendering [=====================>----------------------] at 3.2 fps ~ eta: 16s
## Rendering [=====================>----------------------] at 3.2 fps ~ eta: 15s
## Rendering [======================>---------------------] at 3.2 fps ~ eta: 15s
## Rendering [======================>---------------------] at 3.3 fps ~ eta: 14s
## Rendering [=======================>--------------------] at 3.3 fps ~ eta: 14s
## Rendering [========================>-------------------] at 3.3 fps ~ eta: 13s
## Rendering [=========================>------------------] at 3.3 fps ~ eta: 13s
## Rendering [=========================>------------------] at 3.3 fps ~ eta: 12s
## Rendering [=========================>------------------] at 3.4 fps ~ eta: 12s
## Rendering [==========================>-----------------] at 3.4 fps ~ eta: 12s
## Rendering [==========================>-----------------] at 3.4 fps ~ eta: 11s
## Rendering [===========================>----------------] at 3.4 fps ~ eta: 11s
## Rendering [===========================>----------------] at 3.4 fps ~ eta: 10s
## Rendering [============================>---------------] at 3.5 fps ~ eta: 10s
## Rendering [============================>---------------] at 3.5 fps ~ eta: 9s
## Rendering [=============================>--------------] at 3.5 fps ~ eta: 9s
## Rendering [==============================>-------------] at 3.5 fps ~ eta: 8s
## Rendering [==============================>-------------] at 3.6 fps ~ eta: 8s
## Rendering [===============================>------------] at 3.6 fps ~ eta: 8s
## Rendering [================================>-----------] at 3.6 fps ~ eta: 7s
## Rendering [=================================>----------] at 3.7 fps ~ eta: 6s
## Rendering [==================================>---------] at 3.7 fps ~ eta: 6s
## Rendering [==================================>---------] at 3.7 fps ~ eta: 5s
## Rendering [===================================>--------] at 3.7 fps ~ eta: 5s
## Rendering [====================================>-------] at 3.7 fps ~ eta: 5s
## Rendering [====================================>-------] at 3.7 fps ~ eta: 4s
## Rendering [=====================================>------] at 3.8 fps ~ eta: 4s
## Rendering [=====================================>------] at 3.8 fps ~ eta: 3s
## Rendering [======================================>-----] at 3.8 fps ~ eta: 3s
## Rendering [=======================================>----] at 3.8 fps ~ eta: 3s
## Rendering [=======================================>----] at 3.8 fps ~ eta: 2s
## Rendering [========================================>---] at 3.9 fps ~ eta: 2s
## Rendering [=========================================>--] at 3.9 fps ~ eta: 1s
## Rendering [==========================================>-] at 3.9 fps ~ eta: 1s
## Rendering [===========================================>] at 3.9 fps ~ eta: 0s
## Rendering [============================================] at 3.9 fps ~ eta: 0s

show 1000 draws after the warm-up

labs2 <- c('Draws', '90% HPD')
ggplot() +
  geom_point(data = dfs[1:1000,],
             aes(th1, th2, color = '1'), alpha = 0.3) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '2'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('steelblue', 'blue'), labels = labs2) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA), linetype = c(0, 1), alpha = c(1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

show 4500 draws after the warm-up

labs2 <- c('Draws', '90% HPD')
ggplot() +
  geom_point(data = dfs,
             aes(th1, th2, color = '1'), alpha = 0.3) +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '2'), level = 0.9) +
  coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
  labs(x = 'theta1', y = 'theta2') +
  scale_color_manual(values = c('steelblue', 'blue'), labels = labs2) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA), linetype = c(0, 1), alpha = c(1, 1)))) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Convergence diagnostics

samp <- tt
dim(samp) <- c(dim(tt),1)
samp <- aperm(samp, c(1, 3, 2))
res<-monitor(samp, probs = c(0.25, 0.5, 0.75), digits_summary = 2)
## Inference for the input samples (1 chains: each with iter = 5000; warmup = 2500):
## 
##      Q5  Q50 Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## V1 -1.8 -0.5 0.8 -0.5 0.8  1.00       45       34
## V2 -1.9 -0.4 0.8 -0.5 0.8  1.01       45       33
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).
neff <- res[,'n_eff']
# both theta have owen neff, but for plotting these are so close to each
# other, so that single relative efficiency value is used
reff <- mean(neff/(s/2))

Visual convergence diagnostics

Collapse the data frame with row numbers augmented into key-value pairs for visualizing the chains

dfb <- dfs
sb <- s-warm
dfch <- within(dfb, iter <- 1:sb) %>% gather(grp, value, -iter)

Another data frame for visualizing the estimate of the autocorrelation function

nlags <- 100
dfa <- sapply(dfb, function(x) acf(x, lag.max = nlags, plot = F)$acf) %>%
  data.frame(iter = 0:(nlags)) %>% gather(grp, value, -iter)

A third data frame to visualize the cumulative averages and the 95% intervals

dfca <- (cumsum(dfb) / (1:sb)) %>%
  within({iter <- 1:sb
          uppi <-  1.96/sqrt(1:sb)
          upp <- 1.96/(sqrt(1:sb*reff))}) %>%
  gather(grp, value, -iter)

Visualize the chains

ggplot(data = dfch) +
  geom_line(aes(iter, value, color = grp)) +
  labs(title = 'Trends') +
  scale_color_discrete(labels = c('theta1','theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the autocorrelation function

ggplot(data = dfa) +
  geom_line(aes(iter, value, color = grp)) +
  geom_hline(aes(yintercept = 0)) +
  labs(title = 'Autocorrelation function') +
  scale_color_discrete(labels = c('theta1', 'theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the Monte Carlo error estimates

# labels
labs3 <- c('theta1', 'theta2',
           '95% interval for MCMC error',
           '95% interval for independent MC')
ggplot() +
  geom_line(data = dfca, aes(iter, value, color = grp, linetype = grp)) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb*reff)), linetype = 2) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb)), linetype = 3) +
  geom_hline(aes(yintercept = 0)) +
  coord_cartesian(ylim = c(-1.5, 1.5), xlim = c(0,4000)) +
  labs(title = 'Cumulative averages') +
  scale_color_manual(values = c('red','blue',rep('black', 2)), labels = labs3) +
  scale_linetype_manual(values = c(1, 1, 2, 3), labels = labs3) +
  theme(legend.position = 'bottom', legend.title = element_blank())

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDExLjIiCmF1dGhvcjogIkFraSBWZWh0YXJpLCBNYXJrdXMgUGFhc2luaWVtaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQojIyBNZXRyb3BvbGlzIGFsZ29yaXRobQoKZ2dwbG90MiBpcyB1c2VkIGZvciBwbG90dGluZywgdGlkeXIgZm9yIG1hbmlwdWxhdGluZyBkYXRhIGZyYW1lcwoKYGBge3Igc2V0dXAsIG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGdncGxvdDIpCnRoZW1lX3NldCh0aGVtZV9taW5pbWFsKCkpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkoZ2dhbmltYXRlKQpsaWJyYXJ5KGdnZm9yY2UpCmxpYnJhcnkoTUFTUykKbGlicmFyeShycHJvanJvb3QpCmxpYnJhcnkocnN0YW4pCnJvb3Q8LWhhc19maWxlKCIuQkRBX1JfZGVtb3Nfcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpgYGAKClBhcmFtZXRlcnMgb2YgYSBub3JtYWwgZGlzdHJpYnV0aW9uIHVzZWQgYXMgYSB0b3kgdGFyZ2V0IGRpc3RyaWJ1dGlvbgoKYGBge3IgfQp5MSA8LSAwCnkyIDwtIDAKciA8LSAwLjgKUyA8LSBkaWFnKDIpClNbMSwgMl0gPC0gcgpTWzIsIDFdIDwtIHIKYGBgCgpNZXRyb3BvbGlzIHByb3Bvc2FsIGRpc3RyaWJ1dGlvbiBzY2FsZQoKYGBge3IgfQpzcCA8LSAwLjMKYGBgCgpTYW1wbGUgZnJvbSB0aGUgdG95IGRpc3RyaWJ1dGlvbiB0byB2aXN1YWxpemUgOTAlIEhQRAppbnRlcnZhbCB3aXRoIGdncGxvdCdzIHN0YXRfZWxsaXBzZSgpCgpgYGB7ciB9CmRmdCA8LSBkYXRhLmZyYW1lKG12cm5vcm0oMTAwMDAwLCBjKDAsIDApLCBTKSkKYGBgCgpzZWUgQkRBMyBwLiA4NSBmb3IgaG93IHRvIGNvbXB1dGUgSFBEIGZvciBtdWx0aXZhcmlhdGUgbm9ybWFsCmluIDJkLWNhc2UgY29udG91ciBmb3IgOTAlIEhQRCBpcyBhbiBlbGxpcHNlLCB3aG9zZSBzZW1pbWFqb3IKYXhlcyBjYW4gYmUgY29tcHV0ZWQgZnJvbSB0aGUgZWlnZW52YWx1ZXMgb2YgdGhlIGNvdmFyaWFuY2UKbWF0cml4IHNjYWxlZCBieSBhIHZhbHVlIHNlbGVjdGVkIHRvIGdldCBlbGxpcHNlIG1hdGNoIHRoZQpkZW5zaXR5IGF0IHRoZSBlZGdlIG9mIDkwJSBIUEQuIEFuZ2xlIG9mIHRoZSBlbGxpcHNlIGNvdWxkIGJlCmNvbXB1dGVkIGZyb20gdGhlIGVpZ2VudmVjdG9ycywgYnV0IHNpbmNlIHRoZSBtYXJnaW5hbHMgYXJlIHNhbWUKd2Uga25vdyB0aGF0IGFuZ2xlIGlzIHBpLzQKU3RhcnRpbmcgdmFsdWUgb2YgdGhlIGNoYWluCgpgYGB7ciB9CnQxIDwtIC0yLjUKdDIgPC0gMi41CmBgYAoKTnVtYmVyIG9mIGl0ZXJhdGlvbnMuCgpgYGB7ciB9Ck0gPC0gNTAwMApgYGAKCkluc2VydCB5b3VyIG93biBNZXRyb3BvbGlzIHNhbXBsaW5nIGhlcmUKCmBgYHtyIH0KIyBBbGxvY2F0ZSBtZW1vcnkgZm9yIHRoZSBzYW1wbGUKdHQgPC0gbWF0cml4KHJlcCgwLCAyKk0pLCBuY29sID0gMikKdHRbMSxdIDwtIGModDEsIHQyKSAgICAjIFNhdmUgc3RhcnRpbmcgcG9pbnQKIyBGb3IgZGVtb25zdHJhdGlvbiBsb2FkIHByZS1jb21wdXRlZCB2YWx1ZXMKIyBSZXBsYWNlIHRoaXMgd2l0aCB5b3VyIGFsZ29yaXRobSEKIyB0dCBpcyBhIE0geCAyIGFycmF5LCB3aXRoIE0gZHJhd3Mgb2YgYm90aCB0aGV0YV8xIGFuZCB0aGV0YV8yCmxvYWQocm9vdCgiZGVtb3NfY2gxMSIsImRlbW8xMV8yYS5SRGF0YSIpKQpgYGAKClRoZSByZXN0IGlzIGZvciBpbGx1c3RyYXRpb24KVGFrZSB0aGUgZmlyc3QgMjAwIGRyYXdzCnRvIGlsbHVzdHJhdGUgaG93IHRoZSBzYW1wbGVyIHdvcmtzCgpgYGB7ciB9CmRmMTAwIDwtIGRhdGEuZnJhbWUoaWQ9cmVwKDEsMTAwKSwKICAgICAgICAgICAgICAgICAgICBpdGVyPTE6MTAwLCAKICAgICAgICAgICAgICAgICAgICB0aDEgPSB0dFsxOjEwMCwgMV0sCiAgICAgICAgICAgICAgICAgICAgdGgyID0gdHRbMToxMDAsIDJdLAogICAgICAgICAgICAgICAgICAgIHRoMWwgPSBjKHR0WzEsIDFdLCB0dFsxOigxMDAtMSksIDFdKSwKICAgICAgICAgICAgICAgICAgICB0aDJsID0gYyh0dFsxLCAyXSwgdHRbMTooMTAwLTEpLCAyXSkpCmBgYAoKVGFrZSB0aGUgZmlyc3QgNTAwMCBvYnNlcnZhdGlvbnMgYWZ0ZXIgd2FybXVwIG9mIDUwCgpgYGB7ciB9CnMgPC0gNTAwMAp3YXJtIDwtIDUwMApkZnMgPC0gZGF0YS5mcmFtZSh0aDEgPSB0dFsod2FybSsxKTpzLCAxXSwgdGgyID0gdHRbKHdhcm0rMSk6cywgMl0pCmBgYAoKUmVtb3ZlIHdhcm0tdXAgcGVyaW9kIG9mIDUwIGZpcnN0IGRyYXdzIGxhdGVyCgpgYGB7ciB9CiMgbGFiZWxzIGFuZCBmcmFtZSBpbmRpY2VzIGZvciB0aGUgcGxvdApsYWJzMSA8LSBjKCdEcmF3cycsICdTdGVwcyBvZiB0aGUgc2FtcGxlcicsICc5MCUgSFBEJykKcDEgPC0gZ2dwbG90KCkgKwogIGdlb21faml0dGVyKGRhdGEgPSBkZjEwMCwgd2lkdGg9MC4wNSwgaGVpZ2h0PTAuMDUsCiAgICAgICAgICAgICAgYWVzKHRoMSwgdGgyLCBncm91cD1pZCwgY29sb3IgPScxJyksIGFscGhhPTAuMykgKwogIGdlb21fc2VnbWVudChkYXRhID0gZGYxMDAsIGFlcyh4ID0gdGgxLCB4ZW5kID0gdGgxbCwgY29sb3IgPSAnMicsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHkgPSB0aDIsIHllbmQgPSB0aDJsKSkgKwogIHN0YXRfZWxsaXBzZShkYXRhID0gZGZ0LCBhZXMoeCA9IFgxLCB5ID0gWDIsIGNvbG9yID0gJzMnKSwgbGV2ZWwgPSAwLjkpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoLTQsIDQpLCB5bGltID0gYygtNCwgNCkpICsKICBsYWJzKHggPSAndGhldGExJywgeSA9ICd0aGV0YTInKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoJ3JlZCcsICdmb3Jlc3RncmVlbicsJ2JsdWUnKSwgbGFiZWxzID0gbGFiczEpICsKICBndWlkZXMoY29sb3IgPSBndWlkZV9sZWdlbmQob3ZlcnJpZGUuYWVzID0gbGlzdCgKICAgIHNoYXBlID0gYygxNiwgTkEsIE5BKSwgbGluZXR5cGUgPSBjKDAsIDEsIDEpKSkpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKClRoZSBmb2xsb3dpbmcgZ2VuZXJhdGVzIGEgZ2lmIGFuaW1hdGlvbgpvZiB0aGUgc3RlcHMgb2YgdGhlIHNhbXBsZXIgKG1pZ2h0IHRha2UgMTAgc2Vjb25kcykuCgpgYGB7ciBNZXRyb3BvbGlzICgxKX0KYW5pbWF0ZShwMSArICAgCiAgICAgICAgICB0cmFuc2l0aW9uX3JldmVhbChhbG9uZz1pdGVyKSArIAogICAgICAgICAgc2hhZG93X3RyYWlsKDAuMDEpKQpgYGAKClBsb3QgdGhlIGZpbmFsIGZyYW1lCgpgYGB7ciB9CnAxCmBgYAoKc2hvdyAxMDAwIGRyYXdzIGFmdGVyIHRoZSB3YXJtLXVwCgpgYGB7ciB9CmxhYnMyIDwtIGMoJ0RyYXdzJywgJzkwJSBIUEQnKQpnZ3Bsb3QoKSArCiAgZ2VvbV9wb2ludChkYXRhID0gZGZzWzE6MTAwMCxdLAogICAgICAgICAgICAgYWVzKHRoMSwgdGgyLCBjb2xvciA9ICcxJyksIGFscGhhID0gMC4zKSArCiAgc3RhdF9lbGxpcHNlKGRhdGEgPSBkZnQsIGFlcyh4ID0gWDEsIHkgPSBYMiwgY29sb3IgPSAnMicpLCBsZXZlbCA9IDAuOSkgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygtNCwgNCksIHlsaW0gPSBjKC00LCA0KSkgKwogIGxhYnMoeCA9ICd0aGV0YTEnLCB5ID0gJ3RoZXRhMicpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygnc3RlZWxibHVlJywgJ2JsdWUnKSwgbGFiZWxzID0gbGFiczIpICsKICBndWlkZXMoY29sb3IgPSBndWlkZV9sZWdlbmQob3ZlcnJpZGUuYWVzID0gbGlzdCgKICAgIHNoYXBlID0gYygxNiwgTkEpLCBsaW5ldHlwZSA9IGMoMCwgMSksIGFscGhhID0gYygxLCAxKSkpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpzaG93IDQ1MDAgZHJhd3MgYWZ0ZXIgdGhlIHdhcm0tdXAKCmBgYHtyIH0KbGFiczIgPC0gYygnRHJhd3MnLCAnOTAlIEhQRCcpCmdncGxvdCgpICsKICBnZW9tX3BvaW50KGRhdGEgPSBkZnMsCiAgICAgICAgICAgICBhZXModGgxLCB0aDIsIGNvbG9yID0gJzEnKSwgYWxwaGEgPSAwLjMpICsKICBzdGF0X2VsbGlwc2UoZGF0YSA9IGRmdCwgYWVzKHggPSBYMSwgeSA9IFgyLCBjb2xvciA9ICcyJyksIGxldmVsID0gMC45KSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC00LCA0KSwgeWxpbSA9IGMoLTQsIDQpKSArCiAgbGFicyh4ID0gJ3RoZXRhMScsIHkgPSAndGhldGEyJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdzdGVlbGJsdWUnLCAnYmx1ZScpLCBsYWJlbHMgPSBsYWJzMikgKwogIGd1aWRlcyhjb2xvciA9IGd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXMgPSBsaXN0KAogICAgc2hhcGUgPSBjKDE2LCBOQSksIGxpbmV0eXBlID0gYygwLCAxKSwgYWxwaGEgPSBjKDEsIDEpKSkpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCiMjIyBDb252ZXJnZW5jZSBkaWFnbm9zdGljcwoKYGBge3IgfQpzYW1wIDwtIHR0CmRpbShzYW1wKSA8LSBjKGRpbSh0dCksMSkKc2FtcCA8LSBhcGVybShzYW1wLCBjKDEsIDMsIDIpKQpyZXM8LW1vbml0b3Ioc2FtcCwgcHJvYnMgPSBjKDAuMjUsIDAuNSwgMC43NSksIGRpZ2l0c19zdW1tYXJ5ID0gMikKbmVmZiA8LSByZXNbLCduX2VmZiddCiMgYm90aCB0aGV0YSBoYXZlIG93ZW4gbmVmZiwgYnV0IGZvciBwbG90dGluZyB0aGVzZSBhcmUgc28gY2xvc2UgdG8gZWFjaAojIG90aGVyLCBzbyB0aGF0IHNpbmdsZSByZWxhdGl2ZSBlZmZpY2llbmN5IHZhbHVlIGlzIHVzZWQKcmVmZiA8LSBtZWFuKG5lZmYvKHMvMikpCmBgYAoKIyMjIFZpc3VhbCBjb252ZXJnZW5jZSBkaWFnbm9zdGljcwpDb2xsYXBzZSB0aGUgZGF0YSBmcmFtZSB3aXRoIHJvdyBudW1iZXJzIGF1Z21lbnRlZAppbnRvIGtleS12YWx1ZSBwYWlycyBmb3IgdmlzdWFsaXppbmcgdGhlIGNoYWlucwoKYGBge3IgfQpkZmIgPC0gZGZzCnNiIDwtIHMtd2FybQpkZmNoIDwtIHdpdGhpbihkZmIsIGl0ZXIgPC0gMTpzYikgJT4lIGdhdGhlcihncnAsIHZhbHVlLCAtaXRlcikKYGBgCgpBbm90aGVyIGRhdGEgZnJhbWUgZm9yIHZpc3VhbGl6aW5nIHRoZSBlc3RpbWF0ZSBvZgp0aGUgYXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uCgpgYGB7ciB9Cm5sYWdzIDwtIDUwCmRmYSA8LSBzYXBwbHkoZGZiLCBmdW5jdGlvbih4KSBhY2YoeCwgbGFnLm1heCA9IG5sYWdzLCBwbG90ID0gRikkYWNmKSAlPiUKICBkYXRhLmZyYW1lKGl0ZXIgPSAwOihubGFncykpICU+JSBnYXRoZXIoZ3JwLCB2YWx1ZSwgLWl0ZXIpCmBgYAoKQSB0aGlyZCBkYXRhIGZyYW1lIHRvIHZpc3VhbGl6ZSB0aGUgY3VtdWxhdGl2ZSBhdmVyYWdlcwphbmQgdGhlIDk1JSBpbnRlcnZhbHMKCmBgYHtyIH0KZGZjYSA8LSAoY3Vtc3VtKGRmYikgLyAoMTpzYikpICU+JQogIHdpdGhpbih7aXRlciA8LSAxOnNiCiAgdXBwaSA8LSAgMS45Ni9zcXJ0KDE6c2IpCiAgdXBwIDwtIDEuOTYvKHNxcnQoMTpzYipyZWZmKSl9KSAlPiUKICBnYXRoZXIoZ3JwLCB2YWx1ZSwgLWl0ZXIpCmBgYAoKVmlzdWFsaXplIHRoZSBjaGFpbnMKCmBgYHtyIH0KZ2dwbG90KGRhdGEgPSBkZmNoKSArCiAgZ2VvbV9saW5lKGFlcyhpdGVyLCB2YWx1ZSwgY29sb3IgPSBncnApKSArCiAgbGFicyh0aXRsZSA9ICdUcmVuZHMnKSArCiAgc2NhbGVfY29sb3JfZGlzY3JldGUobGFiZWxzID0gYygndGhldGExJywndGhldGEyJykpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKClZpc3VhbGl6ZSB0aGUgZXN0aW1hdGUgb2YgdGhlIGF1dG9jb3JyZWxhdGlvbiBmdW5jdGlvbgoKYGBge3IgfQpnZ3Bsb3QoZGF0YSA9IGRmYSkgKwogIGdlb21fbGluZShhZXMoaXRlciwgdmFsdWUsIGNvbG9yID0gZ3JwKSkgKwogIGdlb21faGxpbmUoYWVzKHlpbnRlcmNlcHQgPSAwKSkgKwogIGxhYnModGl0bGUgPSAnQXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uJykgKwogIHNjYWxlX2NvbG9yX2Rpc2NyZXRlKGxhYmVscyA9IGMoJ3RoZXRhMScsICd0aGV0YTInKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKVmlzdWFsaXplIHRoZSBlc3RpbWF0ZSBvZiB0aGUgTW9udGUgQ2FybG8gZXJyb3IgZXN0aW1hdGVzCgpgYGB7ciB9CiMgbGFiZWxzCmxhYnMzIDwtIGMoJ3RoZXRhMScsICd0aGV0YTInLAogICAgICAgICAgICc5NSUgaW50ZXJ2YWwgZm9yIE1DTUMgZXJyb3InLAogICAgICAgICAgICc5NSUgaW50ZXJ2YWwgZm9yIGluZGVwZW5kZW50IE1DJykKZ2dwbG90KCkgKwogIGdlb21fbGluZShkYXRhID0gZGZjYSwgYWVzKGl0ZXIsIHZhbHVlLCBjb2xvciA9IGdycCwgbGluZXR5cGUgPSBncnApKSArCiAgZ2VvbV9saW5lKGFlcygxOnNiLCAtMS45Ni9zcXJ0KDE6c2IqcmVmZikpLCBsaW5ldHlwZSA9IDIpICsKICBnZW9tX2xpbmUoYWVzKDE6c2IsIC0xLjk2L3NxcnQoMTpzYikpLCBsaW5ldHlwZSA9IDMpICsKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gMCkpICsKICBjb29yZF9jYXJ0ZXNpYW4oeWxpbSA9IGMoLTEuNSwgMS41KSwgeGxpbSA9IGMoMCw0MDAwKSkgKwogIGxhYnModGl0bGUgPSAnQ3VtdWxhdGl2ZSBhdmVyYWdlcycpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygncmVkJywnYmx1ZScscmVwKCdibGFjaycsIDIpKSwgbGFiZWxzID0gbGFiczMpICsKICBzY2FsZV9saW5ldHlwZV9tYW51YWwodmFsdWVzID0gYygxLCAxLCAyLCAzKSwgbGFiZWxzID0gbGFiczMpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKClNhbWUgYWdhaW4gd2l0aCByPTAuOTkKUGFyYW1ldGVycyBvZiBhIG5vcm1hbCBkaXN0cmlidXRpb24gdXNlZCBhcyBhIHRveSB0YXJnZXQgZGlzdHJpYnV0aW9uCgpgYGB7ciB9CnkxIDwtIDAKeTIgPC0gMApyIDwtIDAuOTkKUyA8LSBkaWFnKDIpClNbMSwgMl0gPC0gcgpTWzIsIDFdIDwtIHIKYGBgCgpNZXRyb3BvbGlzIHByb3Bvc2FsIGRpc3RyaWJ1dGlvbiBzY2FsZQoKYGBge3IgfQpzcCA8LSAwLjMKYGBgCgpTYW1wbGUgZnJvbSB0aGUgdG95IGRpc3RyaWJ1dGlvbiB0byB2aXN1YWxpemUgOTAlIEhQRAppbnRlcnZhbCB3aXRoIGdncGxvdCdzIHN0YXRfZWxsaXBzZSgpCgpgYGB7ciB9CmRmdCA8LSBkYXRhLmZyYW1lKG12cm5vcm0oMTAwMDAwLCBjKDAsIDApLCBTKSkKYGBgCgpzZWUgQkRBMyBwLiA4NSBmb3IgaG93IHRvIGNvbXB1dGUgSFBEIGZvciBtdWx0aXZhcmlhdGUgbm9ybWFsCmluIDJkLWNhc2UgY29udG91ciBmb3IgOTAlIEhQRCBpcyBhbiBlbGxpcHNlLCB3aG9zZSBzZW1pbWFqb3IKYXhlcyBjYW4gYmUgY29tcHV0ZWQgZnJvbSB0aGUgZWlnZW52YWx1ZXMgb2YgdGhlIGNvdmFyaWFuY2UKbWF0cml4IHNjYWxlZCBieSBhIHZhbHVlIHNlbGVjdGVkIHRvIGdldCBlbGxpcHNlIG1hdGNoIHRoZQpkZW5zaXR5IGF0IHRoZSBlZGdlIG9mIDkwJSBIUEQuIEFuZ2xlIG9mIHRoZSBlbGxpcHNlIGNvdWxkIGJlCmNvbXB1dGVkIGZyb20gdGhlIGVpZ2VudmVjdG9ycywgYnV0IHNpbmNlIHRoZSBtYXJnaW5hbHMgYXJlIHNhbWUKd2Uga25vdyB0aGF0IGFuZ2xlIGlzIHBpLzQKU3RhcnRpbmcgdmFsdWUgb2YgdGhlIGNoYWluCgpgYGB7ciB9CnQxIDwtIC0yLjUKdDIgPC0gMi41CmBgYAoKTnVtYmVyIG9mIGl0ZXJhdGlvbnMuCgpgYGB7ciB9Ck0gPC0gNTAwMApgYGAKCkluc2VydCB5b3VyIG93biBNZXRyb3BvbGlzIHNhbXBsaW5nIGhlcmUKCmBgYHtyIH0KIyBBbGxvY2F0ZSBtZW1vcnkgZm9yIHRoZSBzYW1wbGUKdHQgPC0gbWF0cml4KHJlcCgwLCAyKk0pLCBuY29sID0gMikKdHRbMSxdIDwtIGModDEsIHQyKSAgICAjIFNhdmUgc3RhcnRpbmcgcG9pbnQKIyBGb3IgZGVtb25zdHJhdGlvbiBsb2FkIHByZS1jb21wdXRlZCB2YWx1ZXMKIyBSZXBsYWNlIHRoaXMgd2l0aCB5b3VyIGFsZ29yaXRobSEKIyB0dCBpcyBhIE0geCAyIGFycmF5LCB3aXRoIE0gZHJhd3Mgb2YgYm90aCB0aGV0YV8xIGFuZCB0aGV0YV8yCmxvYWQocm9vdCgiZGVtb3NfY2gxMSIsImRlbW8xMV8yYi5SRGF0YSIpKQpgYGAKClRoZSByZXN0IGlzIGZvciBpbGx1c3RyYXRpb24KVGFrZSB0aGUgZmlyc3QgMjAwIGRyYXdzCnRvIGlsbHVzdHJhdGUgaG93IHRoZSBzYW1wbGVyIHdvcmtzCgpgYGB7ciB9CmRmMTAwIDwtIGRhdGEuZnJhbWUoaWQ9cmVwKDEsMTAwKSwKICAgICAgICAgICAgICAgICAgICBpdGVyPTE6MTAwLCAKICAgICAgICAgICAgICAgICAgICB0aDEgPSB0dFsxOjEwMCwgMV0sCiAgICAgICAgICAgICAgICAgICAgdGgyID0gdHRbMToxMDAsIDJdLAogICAgICAgICAgICAgICAgICAgIHRoMWwgPSBjKHR0WzEsIDFdLCB0dFsxOigxMDAtMSksIDFdKSwKICAgICAgICAgICAgICAgICAgICB0aDJsID0gYyh0dFsxLCAyXSwgdHRbMTooMTAwLTEpLCAyXSkpCmBgYAoKVGFrZSB0aGUgZmlyc3QgNTAwMCBvYnNlcnZhdGlvbnMgYWZ0ZXIgd2FybXVwIG9mIDUwCgpgYGB7ciB9CnMgPC0gNTAwMAp3YXJtIDwtIDUwMApkZnMgPC0gZGF0YS5mcmFtZSh0aDEgPSB0dFsod2FybSsxKTpzLCAxXSwgdGgyID0gdHRbKHdhcm0rMSk6cywgMl0pCmBgYAoKUmVtb3ZlIHdhcm0tdXAgcGVyaW9kIG9mIDUwIGZpcnN0IGRyYXdzIGxhdGVyCgpgYGB7ciB9CiMgbGFiZWxzIGFuZCBmcmFtZSBpbmRpY2VzIGZvciB0aGUgcGxvdApsYWJzMSA8LSBjKCdEcmF3cycsICdTdGVwcyBvZiB0aGUgc2FtcGxlcicsICc5MCUgSFBEJykKcDEgPC0gZ2dwbG90KCkgKwogIGdlb21faml0dGVyKGRhdGEgPSBkZjEwMCwgd2lkdGg9MC4wNSwgaGVpZ2h0PTAuMDUsCiAgICAgICAgICAgICBhZXModGgxLCB0aDIsIGdyb3VwPWlkLCBjb2xvciA9JzEnKSwgYWxwaGE9MC4zKSArCiAgZ2VvbV9zZWdtZW50KGRhdGEgPSBkZjEwMCwgYWVzKHggPSB0aDEsIHhlbmQgPSB0aDFsLCBjb2xvciA9ICcyJywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeSA9IHRoMiwgeWVuZCA9IHRoMmwpKSArCiAgc3RhdF9lbGxpcHNlKGRhdGEgPSBkZnQsIGFlcyh4ID0gWDEsIHkgPSBYMiwgY29sb3IgPSAnMycpLCBsZXZlbCA9IDAuOSkgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygtNCwgNCksIHlsaW0gPSBjKC00LCA0KSkgKwogIGxhYnMoeCA9ICd0aGV0YTEnLCB5ID0gJ3RoZXRhMicpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygncmVkJywgJ2ZvcmVzdGdyZWVuJywnYmx1ZScpLCBsYWJlbHMgPSBsYWJzMSkgKwogIGd1aWRlcyhjb2xvciA9IGd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXMgPSBsaXN0KAogICAgc2hhcGUgPSBjKDE2LCBOQSwgTkEpLCBsaW5ldHlwZSA9IGMoMCwgMSwgMSkpKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKVGhlIGZvbGxvd2luZyBnZW5lcmF0ZXMgYSBnaWYgYW5pbWF0aW9uCm9mIHRoZSBzdGVwcyBvZiB0aGUgc2FtcGxlciAobWlnaHQgdGFrZSAxMCBzZWNvbmRzKS4KCmBgYHtyIE1ldHJvcG9saXMgKDIpfQphbmltYXRlKHAxICsgICAKICAgICAgICAgIHRyYW5zaXRpb25fcmV2ZWFsKGFsb25nPWl0ZXIpICsgCiAgICAgICAgc2hhZG93X3RyYWlsKDAuMDEpKQpgYGAKClBsb3QgdGhlIGZpbmFsIGZyYW1lCgpgYGB7ciB9CnAxCmBgYAoKc2hvdyAxMDAwIGRyYXdzIGFmdGVyIHRoZSB3YXJtLXVwCgpgYGB7ciB9CmxhYnMyIDwtIGMoJ0RyYXdzJywgJzkwJSBIUEQnKQpnZ3Bsb3QoKSArCiAgZ2VvbV9wb2ludChkYXRhID0gZGZzWzE6MTAwMCxdLAogICAgICAgICAgICAgYWVzKHRoMSwgdGgyLCBjb2xvciA9ICcxJyksIGFscGhhID0gMC4zKSArCiAgc3RhdF9lbGxpcHNlKGRhdGEgPSBkZnQsIGFlcyh4ID0gWDEsIHkgPSBYMiwgY29sb3IgPSAnMicpLCBsZXZlbCA9IDAuOSkgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygtNCwgNCksIHlsaW0gPSBjKC00LCA0KSkgKwogIGxhYnMoeCA9ICd0aGV0YTEnLCB5ID0gJ3RoZXRhMicpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygnc3RlZWxibHVlJywgJ2JsdWUnKSwgbGFiZWxzID0gbGFiczIpICsKICBndWlkZXMoY29sb3IgPSBndWlkZV9sZWdlbmQob3ZlcnJpZGUuYWVzID0gbGlzdCgKICAgIHNoYXBlID0gYygxNiwgTkEpLCBsaW5ldHlwZSA9IGMoMCwgMSksIGFscGhhID0gYygxLCAxKSkpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpzaG93IDQ1MDAgZHJhd3MgYWZ0ZXIgdGhlIHdhcm0tdXAKCmBgYHtyIH0KbGFiczIgPC0gYygnRHJhd3MnLCAnOTAlIEhQRCcpCmdncGxvdCgpICsKICBnZW9tX3BvaW50KGRhdGEgPSBkZnMsCiAgICAgICAgICAgICBhZXModGgxLCB0aDIsIGNvbG9yID0gJzEnKSwgYWxwaGEgPSAwLjMpICsKICBzdGF0X2VsbGlwc2UoZGF0YSA9IGRmdCwgYWVzKHggPSBYMSwgeSA9IFgyLCBjb2xvciA9ICcyJyksIGxldmVsID0gMC45KSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC00LCA0KSwgeWxpbSA9IGMoLTQsIDQpKSArCiAgbGFicyh4ID0gJ3RoZXRhMScsIHkgPSAndGhldGEyJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdzdGVlbGJsdWUnLCAnYmx1ZScpLCBsYWJlbHMgPSBsYWJzMikgKwogIGd1aWRlcyhjb2xvciA9IGd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXMgPSBsaXN0KAogICAgc2hhcGUgPSBjKDE2LCBOQSksIGxpbmV0eXBlID0gYygwLCAxKSwgYWxwaGEgPSBjKDEsIDEpKSkpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCiMjIyBDb252ZXJnZW5jZSBkaWFnbm9zdGljcwoKYGBge3IgfQpzYW1wIDwtIHR0CmRpbShzYW1wKSA8LSBjKGRpbSh0dCksMSkKc2FtcCA8LSBhcGVybShzYW1wLCBjKDEsIDMsIDIpKQpyZXM8LW1vbml0b3Ioc2FtcCwgcHJvYnMgPSBjKDAuMjUsIDAuNSwgMC43NSksIGRpZ2l0c19zdW1tYXJ5ID0gMikKbmVmZiA8LSByZXNbLCduX2VmZiddCiMgYm90aCB0aGV0YSBoYXZlIG93ZW4gbmVmZiwgYnV0IGZvciBwbG90dGluZyB0aGVzZSBhcmUgc28gY2xvc2UgdG8gZWFjaAojIG90aGVyLCBzbyB0aGF0IHNpbmdsZSByZWxhdGl2ZSBlZmZpY2llbmN5IHZhbHVlIGlzIHVzZWQKcmVmZiA8LSBtZWFuKG5lZmYvKHMvMikpCmBgYAoKIyMjIFZpc3VhbCBjb252ZXJnZW5jZSBkaWFnbm9zdGljcwpDb2xsYXBzZSB0aGUgZGF0YSBmcmFtZSB3aXRoIHJvdyBudW1iZXJzIGF1Z21lbnRlZAppbnRvIGtleS12YWx1ZSBwYWlycyBmb3IgdmlzdWFsaXppbmcgdGhlIGNoYWlucwoKYGBge3IgfQpkZmIgPC0gZGZzCnNiIDwtIHMtd2FybQpkZmNoIDwtIHdpdGhpbihkZmIsIGl0ZXIgPC0gMTpzYikgJT4lIGdhdGhlcihncnAsIHZhbHVlLCAtaXRlcikKYGBgCgpBbm90aGVyIGRhdGEgZnJhbWUgZm9yIHZpc3VhbGl6aW5nIHRoZSBlc3RpbWF0ZSBvZgp0aGUgYXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uCgpgYGB7ciB9Cm5sYWdzIDwtIDEwMApkZmEgPC0gc2FwcGx5KGRmYiwgZnVuY3Rpb24oeCkgYWNmKHgsIGxhZy5tYXggPSBubGFncywgcGxvdCA9IEYpJGFjZikgJT4lCiAgZGF0YS5mcmFtZShpdGVyID0gMDoobmxhZ3MpKSAlPiUgZ2F0aGVyKGdycCwgdmFsdWUsIC1pdGVyKQpgYGAKCkEgdGhpcmQgZGF0YSBmcmFtZSB0byB2aXN1YWxpemUgdGhlIGN1bXVsYXRpdmUgYXZlcmFnZXMKYW5kIHRoZSA5NSUgaW50ZXJ2YWxzCgpgYGB7ciB9CmRmY2EgPC0gKGN1bXN1bShkZmIpIC8gKDE6c2IpKSAlPiUKICB3aXRoaW4oe2l0ZXIgPC0gMTpzYgogICAgICAgICAgdXBwaSA8LSAgMS45Ni9zcXJ0KDE6c2IpCiAgICAgICAgICB1cHAgPC0gMS45Ni8oc3FydCgxOnNiKnJlZmYpKX0pICU+JQogIGdhdGhlcihncnAsIHZhbHVlLCAtaXRlcikKYGBgCgpWaXN1YWxpemUgdGhlIGNoYWlucwoKYGBge3IgfQpnZ3Bsb3QoZGF0YSA9IGRmY2gpICsKICBnZW9tX2xpbmUoYWVzKGl0ZXIsIHZhbHVlLCBjb2xvciA9IGdycCkpICsKICBsYWJzKHRpdGxlID0gJ1RyZW5kcycpICsKICBzY2FsZV9jb2xvcl9kaXNjcmV0ZShsYWJlbHMgPSBjKCd0aGV0YTEnLCd0aGV0YTInKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKVmlzdWFsaXplIHRoZSBlc3RpbWF0ZSBvZiB0aGUgYXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uCgpgYGB7ciB9CmdncGxvdChkYXRhID0gZGZhKSArCiAgZ2VvbV9saW5lKGFlcyhpdGVyLCB2YWx1ZSwgY29sb3IgPSBncnApKSArCiAgZ2VvbV9obGluZShhZXMoeWludGVyY2VwdCA9IDApKSArCiAgbGFicyh0aXRsZSA9ICdBdXRvY29ycmVsYXRpb24gZnVuY3Rpb24nKSArCiAgc2NhbGVfY29sb3JfZGlzY3JldGUobGFiZWxzID0gYygndGhldGExJywgJ3RoZXRhMicpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpWaXN1YWxpemUgdGhlIGVzdGltYXRlIG9mIHRoZSBNb250ZSBDYXJsbyBlcnJvciBlc3RpbWF0ZXMKCmBgYHtyIH0KIyBsYWJlbHMKbGFiczMgPC0gYygndGhldGExJywgJ3RoZXRhMicsCiAgICAgICAgICAgJzk1JSBpbnRlcnZhbCBmb3IgTUNNQyBlcnJvcicsCiAgICAgICAgICAgJzk1JSBpbnRlcnZhbCBmb3IgaW5kZXBlbmRlbnQgTUMnKQpnZ3Bsb3QoKSArCiAgZ2VvbV9saW5lKGRhdGEgPSBkZmNhLCBhZXMoaXRlciwgdmFsdWUsIGNvbG9yID0gZ3JwLCBsaW5ldHlwZSA9IGdycCkpICsKICBnZW9tX2xpbmUoYWVzKDE6c2IsIC0xLjk2L3NxcnQoMTpzYipyZWZmKSksIGxpbmV0eXBlID0gMikgKwogIGdlb21fbGluZShhZXMoMTpzYiwgLTEuOTYvc3FydCgxOnNiKSksIGxpbmV0eXBlID0gMykgKwogIGdlb21faGxpbmUoYWVzKHlpbnRlcmNlcHQgPSAwKSkgKwogIGNvb3JkX2NhcnRlc2lhbih5bGltID0gYygtMS41LCAxLjUpLCB4bGltID0gYygwLDQwMDApKSArCiAgbGFicyh0aXRsZSA9ICdDdW11bGF0aXZlIGF2ZXJhZ2VzJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdyZWQnLCdibHVlJyxyZXAoJ2JsYWNrJywgMikpLCBsYWJlbHMgPSBsYWJzMykgKwogIHNjYWxlX2xpbmV0eXBlX21hbnVhbCh2YWx1ZXMgPSBjKDEsIDEsIDIsIDMpLCBsYWJlbHMgPSBsYWJzMykgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKU2FtZSBhZ2FpbiB3aXRoIHNwID0gMS41CgpgYGB7ciB9CnNwID0gMS41CmBgYAoKSW5zZXJ0IHlvdXIgb3duIE1ldHJvcG9saXMgc2FtcGxpbmcgaGVyZQoKYGBge3IgfQojIEFsbG9jYXRlIG1lbW9yeSBmb3IgdGhlIHNhbXBsZQp0dCA8LSBtYXRyaXgocmVwKDAsIDIqTSksIG5jb2wgPSAyKQp0dFsxLF0gPC0gYyh0MSwgdDIpICAgICMgU2F2ZSBzdGFydGluZyBwb2ludAojIEZvciBkZW1vbnN0cmF0aW9uIGxvYWQgcHJlLWNvbXB1dGVkIHZhbHVlcwojIFJlcGxhY2UgdGhpcyB3aXRoIHlvdXIgYWxnb3JpdGhtIQojIHR0IGlzIGEgTSB4IDIgYXJyYXksIHdpdGggTSBkcmF3cyBvZiBib3RoIHRoZXRhXzEgYW5kIHRoZXRhXzIKbG9hZChyb290KCJkZW1vc19jaDExIiwiZGVtbzExXzJjLlJEYXRhIikpCmBgYAoKVGhlIHJlc3QgaXMgZm9yIGlsbHVzdHJhdGlvbgpUYWtlIHRoZSBmaXJzdCAyMDAgZHJhd3MKdG8gaWxsdXN0cmF0ZSBob3cgdGhlIHNhbXBsZXIgd29ya3MKCmBgYHtyIH0KZGYxMDAgPC0gZGF0YS5mcmFtZShpZD1yZXAoMSwxMDApLAogICAgICAgICAgICAgICAgICAgIGl0ZXI9MToxMDAsIAogICAgICAgICAgICAgICAgICAgIHRoMSA9IHR0WzE6MTAwLCAxXSwKICAgICAgICAgICAgICAgICAgICB0aDIgPSB0dFsxOjEwMCwgMl0sCiAgICAgICAgICAgICAgICAgICAgdGgxbCA9IGModHRbMSwgMV0sIHR0WzE6KDEwMC0xKSwgMV0pLAogICAgICAgICAgICAgICAgICAgIHRoMmwgPSBjKHR0WzEsIDJdLCB0dFsxOigxMDAtMSksIDJdKSkKYGBgCgpUYWtlIHRoZSBmaXJzdCA1MDAwIG9ic2VydmF0aW9ucyBhZnRlciB3YXJtdXAgb2YgNTAKCmBgYHtyIH0KcyA8LSA1MDAwCndhcm0gPC0gNTAwCmRmcyA8LSBkYXRhLmZyYW1lKHRoMSA9IHR0Wyh3YXJtKzEpOnMsIDFdLCB0aDIgPSB0dFsod2FybSsxKTpzLCAyXSkKYGBgCgpSZW1vdmUgd2FybS11cCBwZXJpb2Qgb2YgNTAgZmlyc3QgZHJhd3MgbGF0ZXIKCmBgYHtyIH0KIyBsYWJlbHMgYW5kIGZyYW1lIGluZGljZXMgZm9yIHRoZSBwbG90CmxhYnMxIDwtIGMoJ0RyYXdzJywgJ1N0ZXBzIG9mIHRoZSBzYW1wbGVyJywgJzkwJSBIUEQnKQpwMSA8LSBnZ3Bsb3QoKSArCiAgZ2VvbV9qaXR0ZXIoZGF0YSA9IGRmMTAwLCB3aWR0aD0wLjA1LCBoZWlnaHQ9MC4wNSwKICAgICAgICAgICAgIGFlcyh0aDEsIHRoMiwgZ3JvdXA9aWQsIGNvbG9yID0nMScpLCBhbHBoYT0wLjMpICsKICBnZW9tX3NlZ21lbnQoZGF0YSA9IGRmMTAwLCBhZXMoeCA9IHRoMSwgeGVuZCA9IHRoMWwsIGNvbG9yID0gJzInLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5ID0gdGgyLCB5ZW5kID0gdGgybCkpICsKICBzdGF0X2VsbGlwc2UoZGF0YSA9IGRmdCwgYWVzKHggPSBYMSwgeSA9IFgyLCBjb2xvciA9ICczJyksIGxldmVsID0gMC45KSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC00LCA0KSwgeWxpbSA9IGMoLTQsIDQpKSArCiAgbGFicyh4ID0gJ3RoZXRhMScsIHkgPSAndGhldGEyJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdyZWQnLCAnZm9yZXN0Z3JlZW4nLCdibHVlJyksIGxhYmVscyA9IGxhYnMxKSArCiAgZ3VpZGVzKGNvbG9yID0gZ3VpZGVfbGVnZW5kKG92ZXJyaWRlLmFlcyA9IGxpc3QoCiAgICBzaGFwZSA9IGMoMTYsIE5BLCBOQSksIGxpbmV0eXBlID0gYygwLCAxLCAxKSkpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpUaGUgZm9sbG93aW5nIGdlbmVyYXRlcyBhIGdpZiBhbmltYXRpb24Kb2YgdGhlIHN0ZXBzIG9mIHRoZSBzYW1wbGVyIChtaWdodCB0YWtlIDEwIHNlY29uZHMpLgoKYGBge3IgTWV0cm9wb2xpcyAoMyl9CmFuaW1hdGUocDEgKyAgIAogICAgICAgICAgdHJhbnNpdGlvbl9yZXZlYWwoYWxvbmc9aXRlcikgKyAKICAgICAgICBzaGFkb3dfdHJhaWwoMC4wMSkpCmBgYAoKc2hvdyAxMDAwIGRyYXdzIGFmdGVyIHRoZSB3YXJtLXVwCgpgYGB7ciB9CmxhYnMyIDwtIGMoJ0RyYXdzJywgJzkwJSBIUEQnKQpnZ3Bsb3QoKSArCiAgZ2VvbV9wb2ludChkYXRhID0gZGZzWzE6MTAwMCxdLAogICAgICAgICAgICAgYWVzKHRoMSwgdGgyLCBjb2xvciA9ICcxJyksIGFscGhhID0gMC4zKSArCiAgc3RhdF9lbGxpcHNlKGRhdGEgPSBkZnQsIGFlcyh4ID0gWDEsIHkgPSBYMiwgY29sb3IgPSAnMicpLCBsZXZlbCA9IDAuOSkgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygtNCwgNCksIHlsaW0gPSBjKC00LCA0KSkgKwogIGxhYnMoeCA9ICd0aGV0YTEnLCB5ID0gJ3RoZXRhMicpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygnc3RlZWxibHVlJywgJ2JsdWUnKSwgbGFiZWxzID0gbGFiczIpICsKICBndWlkZXMoY29sb3IgPSBndWlkZV9sZWdlbmQob3ZlcnJpZGUuYWVzID0gbGlzdCgKICAgIHNoYXBlID0gYygxNiwgTkEpLCBsaW5ldHlwZSA9IGMoMCwgMSksIGFscGhhID0gYygxLCAxKSkpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpzaG93IDQ1MDAgZHJhd3MgYWZ0ZXIgdGhlIHdhcm0tdXAKCmBgYHtyIH0KbGFiczIgPC0gYygnRHJhd3MnLCAnOTAlIEhQRCcpCmdncGxvdCgpICsKICBnZW9tX3BvaW50KGRhdGEgPSBkZnMsCiAgICAgICAgICAgICBhZXModGgxLCB0aDIsIGNvbG9yID0gJzEnKSwgYWxwaGEgPSAwLjMpICsKICBzdGF0X2VsbGlwc2UoZGF0YSA9IGRmdCwgYWVzKHggPSBYMSwgeSA9IFgyLCBjb2xvciA9ICcyJyksIGxldmVsID0gMC45KSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC00LCA0KSwgeWxpbSA9IGMoLTQsIDQpKSArCiAgbGFicyh4ID0gJ3RoZXRhMScsIHkgPSAndGhldGEyJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdzdGVlbGJsdWUnLCAnYmx1ZScpLCBsYWJlbHMgPSBsYWJzMikgKwogIGd1aWRlcyhjb2xvciA9IGd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXMgPSBsaXN0KAogICAgc2hhcGUgPSBjKDE2LCBOQSksIGxpbmV0eXBlID0gYygwLCAxKSwgYWxwaGEgPSBjKDEsIDEpKSkpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCiMjIyBDb252ZXJnZW5jZSBkaWFnbm9zdGljcwoKYGBge3IgfQpzYW1wIDwtIHR0CmRpbShzYW1wKSA8LSBjKGRpbSh0dCksMSkKc2FtcCA8LSBhcGVybShzYW1wLCBjKDEsIDMsIDIpKQpyZXM8LW1vbml0b3Ioc2FtcCwgcHJvYnMgPSBjKDAuMjUsIDAuNSwgMC43NSksIGRpZ2l0c19zdW1tYXJ5ID0gMikKbmVmZiA8LSByZXNbLCduX2VmZiddCiMgYm90aCB0aGV0YSBoYXZlIG93ZW4gbmVmZiwgYnV0IGZvciBwbG90dGluZyB0aGVzZSBhcmUgc28gY2xvc2UgdG8gZWFjaAojIG90aGVyLCBzbyB0aGF0IHNpbmdsZSByZWxhdGl2ZSBlZmZpY2llbmN5IHZhbHVlIGlzIHVzZWQKcmVmZiA8LSBtZWFuKG5lZmYvKHMvMikpCmBgYAoKIyMjIFZpc3VhbCBjb252ZXJnZW5jZSBkaWFnbm9zdGljcwpDb2xsYXBzZSB0aGUgZGF0YSBmcmFtZSB3aXRoIHJvdyBudW1iZXJzIGF1Z21lbnRlZAppbnRvIGtleS12YWx1ZSBwYWlycyBmb3IgdmlzdWFsaXppbmcgdGhlIGNoYWlucwoKYGBge3IgfQpkZmIgPC0gZGZzCnNiIDwtIHMtd2FybQpkZmNoIDwtIHdpdGhpbihkZmIsIGl0ZXIgPC0gMTpzYikgJT4lIGdhdGhlcihncnAsIHZhbHVlLCAtaXRlcikKYGBgCgpBbm90aGVyIGRhdGEgZnJhbWUgZm9yIHZpc3VhbGl6aW5nIHRoZSBlc3RpbWF0ZSBvZgp0aGUgYXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uCgpgYGB7ciB9Cm5sYWdzIDwtIDEwMApkZmEgPC0gc2FwcGx5KGRmYiwgZnVuY3Rpb24oeCkgYWNmKHgsIGxhZy5tYXggPSBubGFncywgcGxvdCA9IEYpJGFjZikgJT4lCiAgZGF0YS5mcmFtZShpdGVyID0gMDoobmxhZ3MpKSAlPiUgZ2F0aGVyKGdycCwgdmFsdWUsIC1pdGVyKQpgYGAKCkEgdGhpcmQgZGF0YSBmcmFtZSB0byB2aXN1YWxpemUgdGhlIGN1bXVsYXRpdmUgYXZlcmFnZXMKYW5kIHRoZSA5NSUgaW50ZXJ2YWxzCgpgYGB7ciB9CmRmY2EgPC0gKGN1bXN1bShkZmIpIC8gKDE6c2IpKSAlPiUKICB3aXRoaW4oe2l0ZXIgPC0gMTpzYgogICAgICAgICAgdXBwaSA8LSAgMS45Ni9zcXJ0KDE6c2IpCiAgICAgICAgICB1cHAgPC0gMS45Ni8oc3FydCgxOnNiKnJlZmYpKX0pICU+JQogIGdhdGhlcihncnAsIHZhbHVlLCAtaXRlcikKYGBgCgpWaXN1YWxpemUgdGhlIGNoYWlucwoKYGBge3IgfQpnZ3Bsb3QoZGF0YSA9IGRmY2gpICsKICBnZW9tX2xpbmUoYWVzKGl0ZXIsIHZhbHVlLCBjb2xvciA9IGdycCkpICsKICBsYWJzKHRpdGxlID0gJ1RyZW5kcycpICsKICBzY2FsZV9jb2xvcl9kaXNjcmV0ZShsYWJlbHMgPSBjKCd0aGV0YTEnLCd0aGV0YTInKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKVmlzdWFsaXplIHRoZSBlc3RpbWF0ZSBvZiB0aGUgYXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uCgpgYGB7ciB9CmdncGxvdChkYXRhID0gZGZhKSArCiAgZ2VvbV9saW5lKGFlcyhpdGVyLCB2YWx1ZSwgY29sb3IgPSBncnApKSArCiAgZ2VvbV9obGluZShhZXMoeWludGVyY2VwdCA9IDApKSArCiAgbGFicyh0aXRsZSA9ICdBdXRvY29ycmVsYXRpb24gZnVuY3Rpb24nKSArCiAgc2NhbGVfY29sb3JfZGlzY3JldGUobGFiZWxzID0gYygndGhldGExJywgJ3RoZXRhMicpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpWaXN1YWxpemUgdGhlIGVzdGltYXRlIG9mIHRoZSBNb250ZSBDYXJsbyBlcnJvciBlc3RpbWF0ZXMKCmBgYHtyIH0KIyBsYWJlbHMKbGFiczMgPC0gYygndGhldGExJywgJ3RoZXRhMicsCiAgICAgICAgICAgJzk1JSBpbnRlcnZhbCBmb3IgTUNNQyBlcnJvcicsCiAgICAgICAgICAgJzk1JSBpbnRlcnZhbCBmb3IgaW5kZXBlbmRlbnQgTUMnKQpnZ3Bsb3QoKSArCiAgZ2VvbV9saW5lKGRhdGEgPSBkZmNhLCBhZXMoaXRlciwgdmFsdWUsIGNvbG9yID0gZ3JwLCBsaW5ldHlwZSA9IGdycCkpICsKICBnZW9tX2xpbmUoYWVzKDE6c2IsIC0xLjk2L3NxcnQoMTpzYipyZWZmKSksIGxpbmV0eXBlID0gMikgKwogIGdlb21fbGluZShhZXMoMTpzYiwgLTEuOTYvc3FydCgxOnNiKSksIGxpbmV0eXBlID0gMykgKwogIGdlb21faGxpbmUoYWVzKHlpbnRlcmNlcHQgPSAwKSkgKwogIGNvb3JkX2NhcnRlc2lhbih5bGltID0gYygtMS41LCAxLjUpLCB4bGltID0gYygwLDQwMDApKSArCiAgbGFicyh0aXRsZSA9ICdDdW11bGF0aXZlIGF2ZXJhZ2VzJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdyZWQnLCdibHVlJyxyZXAoJ2JsYWNrJywgMikpLCBsYWJlbHMgPSBsYWJzMykgKwogIHNjYWxlX2xpbmV0eXBlX21hbnVhbCh2YWx1ZXMgPSBjKDEsIDEsIDIsIDMpLCBsYWJlbHMgPSBsYWJzMykgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoK